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Abstract 

We investigate the statistical mechanics of a torsionally constrained 
polymer. The polymer is modeled as a fluctuating rod with bend stiffness 
AksT and twist stiffness CksT. In such a model, thermal bend fluctuations 
couple geometrically to an apphed torque through the relation Lk = Tw+ Wr. 
We explore this coupling and flnd agreement between the predictions of our 
model and recent experimental results on single A-DNA molecules. This 
analysis affords an experimental determination of the microscopic twist 
stiffness (averaged over a helix repeat). Quantitative agreement between 
theory and experiment is obtained using C = 109 nm [i.e. twist rigidity 
CksT =4.5 X 10~^^ergcm). The theory further predicts a thermal reduction 
of the effective twist rigidity induced by bend fluctuations. Finally, we flnd a 
small reflection of molecular chirality in the experimental data and interpret 
it in terms of a twist-stretch coupling of the DNA duplex. 
PACS: 87.15.-v, 87.10.+e, 87.15.By. 
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I Introduction 



In this paper, we investigate the statistical mechanics of a polymer chain with torsional 
rigidity. We model the polymer as an elastic rod subject to thermal fluctuations. Each 
conformation of the chain is statistically weighted according to the energy associated with 
bending and twisting. This is in contrast to conventional polymer models, which account 
only for the energy cost of bending the polymer backbone.^ This neglect of torsional 
energy is often well justified, as many polymers are free to release twist by swiveling about 
the single carbon bonds that constitute their backbone. Even for polymers that cannot 
swivel freely, the twist usually amounts to an uncoupled Gaussian degree of freedom 
that can simply be integrated away. The situation is quite different, however, in the 
presence of a torsional constraint. In this case, the twist is coupled to the conformation 
of the backbone and cannot be eliminated so easily. Such a situation can arise when the 
polymer is hgated into a circle, or when its ends are clamped and a torque is applied at 
one end. The concept of a torsional constraint can also be extended to the dynamics of 
a polymer in a viscous fluid: here viscous damping provides the necessary resistance to 
the stress.0 Whatever the origin of the constraint, it will result in a coupling between 
the twist and the bending modes of the backbone. 

The origin of this coupling lies in White's theorem: Lk = Tw+ Wr^' ^ This formula 
relates a global topological invariant of any pair of closed curves (the Linking number, 
Lk), to the sum of a local strain field (the Twist, Tw) and a global configurational 
integral (the Writhe, Wr). If the linking number is fixed, the polymer will be forced 
to distribute the invariant Lk between the degrees of freedom associated with Tw and 
Wr. From a statistical mechanics point of view, the set of complexions available to the 
system is then restricted. The elastic energy of each allowed complexion reflects the sum 
of a twisting energy and a bending energy associated with the Writhe of the backbone. 
Of course we do not need to consider fixed linking number for torsional rigidity to be 
important: a chemical potential for Lk in the form of an applied torque also couples the 
bend fluctuations to the twist. 

Perhaps the most important examples of twist-storing polymers are biopolymers, 
especially DNA. Unlike many of its hydrocarbon-chain cousins, the monomers of DNA 
are joined by multiple covalent bonds; additional specific pairing interactions between 
bases prevent slippage between the strands. This multiply-bonded structure inhibits the 
unwinding of the DNA helix to release a torsional stress; instead, there is an elastic 
energy cost associated with the deformation. 

Recently it has become possible to perform experiments on single molecules of DNA. 
In a classic experiment. Smith et a/.E anchored one end of a DNA duplex to a solid 
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substrate while the other end was attached to a magnetic bead. The conformations 
of the polymer could then be probed by considering the end-to-end extension of the 
chain as a function of the magnetic force applied to the bead. These experiments, 
and others which stretch DNA molecules using electric fields,^ hydrodynamic flows,^ 
or optical tweezerJi^ were soon analyzed using the "worm-like chain" (WLC) model.E 
Working within this framework, Bustamante, Marko and Siggi a0 and Vologo dskiiHl 
were able to reproduce the experimental force-extension curves for DNA over a wide 
range of forces (from 10~^ pN to 10 pN) with just one fitting parameter, the DNA bend 
persistence length. 

Since the original DNA stretching experiments, significant improvements have been 



made. In particular, a series of elegant experiments" has succeeded in torsionally 

constraining the DNA using swivel- free attachments at both ends. As a result, one can 
now directly explore the interplay between DNA's internal resistance to twisting and the 
conformations of its backbone. 

In this paper, we will explain some of these new results analytically in terms of a 
theory of twist-storing polymers. Our final formula, given in ( ^T]) below, quantitatively 
fits the experimental data of Strick et a/.0 and of Allemand and Croquette with 
just two important fit parameters: the bend stiffness A and twist stiffness C (a more 
precise statement appears below). Our analytical approach rests upon linear elasticity 
and perturbation theory about a straight rod. Thus we do not address the remarkable 



structural transitions induced in DNA by torsional stress," " nor will we systematically 
study the plectonemic transition or other phenomena involving self- avoidance. Marko and 
Siggia have previously studied the effects of thermal fluctuations on plectonemic DNA;II3 
we have chosen instead to work in a regime not afflicted by this theoretical difflculty. We 
will show that our analysis is justified in a well-defined region of parameter space where 
many experimental data points are available (solid symbols in Figure |l|), and from the 
data deduce the fundamental elastic parameters of DNA. 

The main points of our results were announced previously.0 EH Some of the steps 
were independently derived by Bouchiat and MezardS in a different analysis of the 
same experiments. The present paper gives some new analytical results, particularly in 
section V.D| , and applies the analysis to some new experimental data (see Figure |l|). 

In addition to these analytical results, Vologodskii and Marko, and Bouchiat 
and Mezard, have recently performed Monte Carlo simulationsEH EHI to study the 
conformations of DNA under applied tensions and torques appropriate to those in the 
experiments studied here. Marko has also studied the related problem of torsional 
constraints on the overstretching transition.c^' c2 
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Apart from quantitatively reproducing the experimental extension curves with just a 
few fit parameters, our theory also predicts a reduction of the effective twist rigidity of 
a polymer caused by conformational fluctuations. We give the form of a new effective 
twist rigidity Ccsk^T, which is smaller than the microscopic rigidity Ck^^T. This effect, 



anticipated some time ago by Shimada and Yamakawa^^ has a simple explanation: part 
of the excess Link imposed on a solid rod can be moved into the bend deformations of its 
backbone through the coupling associated with the Lie constraint. Our simple formula 
((|^) below) makes this intuition precise for the case of a highly stretched rod. 

It may at first seem that all the relevant physics could be found in the classical 



works of the nineteenth century,c3 but actually one can see at once that classical beam 
theory is qualitatively at odds with the experimental data of Figure |l|: it says that a 
rod under tension will simply twist in response to an applied torque r as long as r is 
small enough. Only when the torque exceeds a critical value will the rod buckle into a 
helical configuration, thus shortening the end-to-end extension. Unlike its macroscopic 
counterpart, however, a microscopic rod is continuously buffeted by thermal fluctuations. 
Because the rod is never straight, its average shape will respond as soon as any torsional 
stress is applied; there is no threshold, as seen in Figure |l|. In sections p^-[V| we will 
create a simple mathematical model embodying this observation and use it to explain 
the data. 

II Experiment 

The statistical mechanical problem of a twist-storing polymer subject to a Lk constraint 



is realized in the experiments of Strick et al)^ li^ and Allemand and Croquette.B:^ In 
these experiments, a segment of double-stranded A-DNA of length L ^ 15.6 yum is held 
at both ends: one end is fixed to a glass plate while the other is attached to a magnetic 
bead. Both ends are bound in such a way as to prevent swiveling of the polymer about 
the point of attachment. By rotating the magnetic bead in an applied magnetic field, 
the experimenters are then able to adjust the excess linking number to any desired, fixed 
value. 

While the direction of the applied field fixes the linking number, a gradient in the same 
field allows the DNA molecules to be put under tension. The experiment is therefore able 
to study the statistical mechanics of the biopolymer in the fixed tension / and linking 
number Lk ensemble. The measured response is then the end-to-end extension Lk) of 
the chain as a function of the applied stress. In contrast, traditional ligation experiments 
control only L and Lk, and Lk/L can take on only rather widely-spaced discrete values. 
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Moreover, the measured quantity is gel mobility, whose relation to backbone conformation 
is not simple. 

Some of the experimental results for forces greater than 0.1 pN are shown in Figure |l|. 
In the figure, the solid lines are our theoretical fit to the solid points. These curves were 
produced by fitting four parameters: the microscopic persistence lengths A, C and twist- 
stretch coupling D (all averaged over a helical repeat), as well as the arclength of the 
polymer L. The bend persistence length A has been determined in a number of earlier 
experiments,^' HH' while L can be determined from only the data points with zero 
excess Link. The fitted values of A and L therefore serve mainly as a check of the theory. 
In our fit we used 69 different points, only some of which are depicted as the solid symbols 
in Figure |l]. The figure also shows open symbols. These points correspond to {f,Lk) 
pairs that lie outside the region where our model, which has no explicit self-avoidance, 
is valid. Due to this neglect of self-avoidance, our phantom chain model will have a 
mathematical pathology associated with configurations that include self- crossings. To 
deal with these difficulties, we will simply require that the chain be pulled hard enough 
that such configurations become statistically negligible. As we will see, "pulling hard 
enough" corresponds to a restriction on the apphed stretching force / and the applied 
torque r (see appendix B). Apart from the restrictions of the phantom chain model, 
there were also omissions of data points for physical reasons. For example, at large 
applied tensions and torques, the DNA molecule undergoes structural transformations. 
In section |VI|, we will discuss our data selection criteria and the fitting procedure more 
fully. 



Ill Physical Model 

Throughout most of this paper we will model DNA as a fiuctuating elastic rod of uniform 
circular cross-section and fixed contour length L. This idealization neglects DNA's helical 
nature: in particular, the length scale associated with the helical pitch of the molecule 
{2it/ujo = 3.6 nm) does not enter as a parameter. The concept of fractional overtwist 
{a = 2TxALk/ Luoq) is therefore meaningless. Nevertheless, we will retain the traditional 
notation to provide a connection to the published experimental data, expressing our 
results in terms of a and noting that a and ujq enter only in the combination au^. In 
the main text we will show that our achiral, isotropic elastic rod model captures the 
main features of Figure At the end of our calculation, in (p|, we will also allow for 
intrinsic stretching and a possible asymmetry between positive and negative a, a chiral 
effect associated with the twist-stretch coupling of a helical rod. 
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In appendix A we will introduce helical pitch effects and show that at modest 
stretching tension they can be summarized in an effective, "coarse-grained" energy (see 
(|l]) below). They also lead to a new phenomenon, chiral entropic elasticity, via the 
twist-bend coupling of DNA.^ This effect is potentially another source of asymmetry 
between over- and undertwisting, but the available data do not at present give detailed 
information about the asymmetry, and so we omit this complication from the main text. 

Accordingly we define an elastic energy functional which describes the bending and 

p 

twisting of an isotropic elastic rod of fixed arclength Lx 
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- ^ /' {di/dsYds, and = ^ r ds. 

Jo Kr,l 2 Jo 
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In these formulas t{s) is the tangent to the rod backbone at the point with arclength s 
from the end. We imagine inscribing permanently a "material frame" embedded in the 
rod; then ^3 is the rate of rotation of this frame about t (see ([Tl|) below; our notation 



mainly follows that of Marko and Siggiali3) . We are free to choose a convenient material 
frame; we choose one which coincides with the fixed lab frame when the molecule is 
unstressed. (In keeping with the remarks above, there is no reason to choose a material 
frame initially rotating relative to the lab at tuo-) A and C are the bend and twist 
"persistence lengths," which are given by the respective elastic constants divided by 
k^T. These parameters are understood to be averaged (or "coarse-grained") over the 
scale of a helical repeat. In appendix A we find the relation between them and a more 
elaborate elasticity theory incorporating the intrinsic helicity of the DNA duplex. 

Equations (|l]) are mathematically identical to the kinetic energy of a symmetric 
spinning top with arclength s playing the role of time. Hence there is a direct analogy 
between the dynamical equations of motion for a top and the equations describing the 



equilibrium for an elastic rod, an observation due to Kirchoff.l^ The main technical 



point of our analysis is the extension of Kirchoff's observation to a mathematical 
correspondence between the thermal fluctuations of an elastic rod and the quantum 
mechanics of a spinning top. 



20, 19,01 



The bend persistence length A which appears in ([^) is a well-known parameter that 
has been measured in several experiments. Among other things, this parameter is known 
to depend on the salt concentration of the surrounding fiuid.0 Wang et al. have measured 
A = 47 nm for DNA in buffer conditions similar to those in the experiments studied 
here.E 
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The value of the twist persistence length C has not been determined as directly 
as A. Cyclization kinetics studies,0 EH EH topoisomer distribution analyse^H' EH 
and fluorescence polarization anisotropy (FPA) experimentdH' 0' HI have provided 
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measurements of this parameter, but these determinations are somewhat indirect and 
the results have been difficult to reconcile with each other .H HE In particular, results 
obtained from straight and circular DNA's using a single technique (FPA) yield different 
values of the twist rigidity: C ~ 50 nm for linear DNA's and C ~ 85 nm for circular 
DNA's.^ This discrepancy may be a consequence of the thermal softening of the torsional 
rigidity predicted by our theory (see (||)). The main goal of the present paper is to 
interpret the single DNA molecule data in Figure in terms of a theory we call "torsional 
directed walks", thereby permitting a new measurement of C. Like the bending rigidity 
A, C may be expected to depend on the buffer solution; the dependence of C should 
however be much weaker than A since twisting does not modify the spatial distribution 
between charges on the backbone to the same degree as bending. 

The rod is subject to a stretching force / and a torsional constraint. It will prove 
simplest to impose the torsional constraint through a fixed applied torque r rather than 
directly through a fixed linking number. Since the molecules we will study are many 
times longer than A or C, we are in the thermodynamic limit, and so we expect the two 
ensembles to give the same physical results. 

The two stresses on the polymer require the introduction of two more terms in the 
polymer's energy functional: 

-/■z = -/ Tt-e^ds, and = _2'Kf ■ Lk. (2) 



k^T Jo k^T 

Here z is the end-to-end extension of the polymer. The tension and torque have been 
expressed in terms of the thermal energy: 

/ = ///cbT, and f = r/kj,T. (3) 

In (0) and throughout this paper, Lk denotes the excess Link, consistent with the remarks 
at the beginning of this section; thus Lk = for the unstressed rod. In general, Lk is 
defined only for closed loops. If we have an open chain with both ends held at fixed 
orientations, as in the experiments under study, then we can draw a fixed, imaginary 
return path completing our chain to a closed loop and let Lk denote the Link of this 
closed loop. Choosing the return path so that Lk = when the rod is straight and 
unstressed then gives in general Lk = Tw + Wr where the terms on the right refer only 
to the open, physical rod. 

Before we include -Etorque in our energy functional, the Link must be more explicitly 
expressed. To get a useful expression, we first note that the Twist is defined as 

Tw=^ n^ds. (4) 
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The Writhe involves only the space curve f{s) swept out by the rod's centerline. In 
general, this number is given by a complicated, non-local formulaS' involving a 
double integral around the closed curve: 

wr=^idsfis-(-^.'^]. ,:j'^- y;> . (5) 

4:71 J J \ ds ds J |r(s)-r(s')P 

However, a result due to Fuller allows us to rewrite this quantity as a single integral 
over a local Writhe density. This simplification is made possible by noting that for 
small variations about some reference curve ro{s), the integrand in (|^) becomes a total 
derivative. Performing one of the integrals then yields a single integral over a local 
quantity.!! Specializing to the case where the reference curve is just the Cz-axis then 
givesCiJ 



Wr=^ f'-^^^^^S^ds. (6) 

Fuller's result holds as long as there is a continuous set of non-self-intersecting curves 
interpolating between the reference curve and the curve in question, such that the 
denominator in (P) never vanishes. We can now combine the terms to get the full energy 
functional for our model of DNA: 

E i?bend , -Stwist x o ~ r 7 /'7^ 

H — ; J ■ z — zvrr ■ Lk. (7) 
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Formulas |^ define the elastic model we will use through the end of 

section [V.D| . Later, in section and appendix A we will consider various elaborations 



of the model and determine that they are relatively unimportant in capturing the main 
features of the data in Figure |I|. 

As noted in the introduction, we expect that thermal fiuctuations will have an 
important effect on the rod's twist degree of freedom. A macroscopic elastic rod under 
tension will sustain a finite amount of applied torsional stress without buckling. Once 
a threshold is reached, however, the stress can be partially relaxed by bending the 
backbone. Linear stability analysis of the energy (|^) shows that this threshold is given 
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by Tci-it = 2yA/. Unlike its macroscopic counterpart, however, a microscopic rod is 
subject to thermal fiuctuations. These fiuctuations prevent the rod from ever being 
straight; as we show below, even infinitesimal torsional stresses will then affect the bend 
fiuctuations. Even though there is no chiral energy term, individual fiuctuations will 
not be inversion symmetric. An applied torsion will push the fiuctuations with the 
corresponding helical sense closer to instability, while suppressing those of the opposite 
helical sense. The end result will be a coupling between the applied torsion and the mean 
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end-to-end extension of the rod proportional to (terms linear in r must drop out since 
the model does not break inversion symmetry). 

Later we will consider the effects of molecular chirality: e.g. in section |V^, we will 

99 d9 d*^ 

include a twist-stretch coupling term D.\=^ [i=l It will turn out that the effect of this 
coupling on the experiment we study is small: this is already apparent in Figure |1| where 
the data points are nearly symmetric about a = 0. Nevertheless, by including the twist- 
stretch coupling, we will be able to determine the parameter D roughly. 

Another way that chirality enters a physical model of DNA is through an anisotropic 
bending term. Any transverse slice through the molecule is easier to bend in one direction 
than in another. Microscopically, this anisotropy has its origin in the shape of the base 
pair plates that make up the rungs on the DNA ladder. Since these plates are longer 
in one direction than the other, bending about the short axis ("tilt") is more difficult 
than bending about the long axis ("roll").0 El' EH In appendix A we consider such an 
anisotropy, as well as the related twist-bend coupling,^ finding that these effects can be 
summarized to good accuracy in an effective coarse-grained model of the form (p. This 
conclusion could have been anticipated since the important fluctuations are on length 



scales around 2ti^ A/ f, and for the forces below 8 pN that we consider, this averages 
over at least several helical repeats. We conclude that the treatment of DNA as an achiral 
rod of elastic material is sufficient to understand how its extension changes under applied 
tension and torque. 

At this point it may be noted that unstressed natural DNA is not a perfect helix; 
its axial symmetry is already broken, even in the absence of thermal fluctuations. In 
particular, it is well known that the unstressed, zero temperature structure of DNA is 



sequence dependent .EH' E§l The effect of this quenched disorder has been studied recently 



by Bensimon, Dohmi, and Mezardil and by one of us.S For simple models of weak 



disorder, the main effect is simply to renormalize the bend persistence length A. In 
the present paper, we neglect explicit inclusion of the quenched disorder associated with 
sequence-dependent effects. Thus our bend rigidity A is the effective value including 
disorder. 

Even though the bend and twist rigidities represent averages over a helix repeat, they 
are still microscopic parameters and therefore reflect only the short-scale behavior. As we 
go to longer length scales, we expect the effective bend and twist rigidities to be modified 
by the geometric coupling implicit in White's formula. In particular, we will find that 
the effective twist rigidity is reduced for small applied tensions: 

\ -1 

Ceff = C I 1 ^ 

4Aa 
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The dependence of Ces on length scale enters through /: as mentioned above, yA/f 
sets the scale of the most important fluctuations in the problem. At small tensions, or 
equivalently at long length scales, C is effectively reduced. Equation describes this 
"softening" of the twist rigidity. The reducing factor is explicitly dependent on k^T, 
indicating that this is a thermal effect. 



IV Group Language 

In the next section we will consider the thermodynamic complexions available to a 
torsionally constrained polymer. To prepare for the task, we must first define convenient 
variables for evaluating the energy functional of the last section on the group of rotations, 
5*0 (3). The bending and twisting deformations that appear in (|1|) as well as the Lagrange 
multiplier terms for extension and Link which appear in will need to be expressed in 
terms of these variables. 

We will use two reference frames related by an element of the rotation group. The 
first of these frames is "space-fixed"; we will take as its basis the orthonormal triad {cj}, 
with i =x,y, or z. A rotation g(s) relates this frame to the "body-fixed" (or "material") 
frame {Ea{s)} with a = 1,2, or 3, where s denotes a point on the rod backbone. As 
mentioned earlier, we will take E^i^s) = i{s) to be the tangent to the rod's centerline, 
and the remaining two vectors to be constant directions when the rod is straight and 
unstressed. The local orientation of the polymer is then given by the 3x3 orthogonal 
matrix gai{s) = Ea{s) -Cj. The matrix g contains only three independent entries. We will 
sometimes find it convenient to represent it in a nonredundant way using Euler angles: 

g(s) = e''-3V(s)g-Lie(s)g-L30(s)_ 

Thus for example i{s) ■ = g3z(s) = cos6'(s). 

The generators of infinitesimal rotations are then matrix operators acting on g. When 
these operators act from the left they are called "body-fixed rotations"; when they act 
from the right they are called "space-fixed rotations" . In either case a convenient basis 
for the generators is 
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and Lg = I -1 I . (10) 




We can then describe the rotation of the material frame as we walk along the rod backbone 
as an infinitesimal body-fixed rotation Q or as a space-fixed rotation Q, where 



Q = gg-^ and Q = g-'g. (11) 
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Here and elsewhere, a dot signifies d/ds. We will also write the projections of the rotation 
rates onto the generators as 

fi„ = (Q,U) = -lTr[QU] (12) 

and similarly for Cli. 

With these definitions we can cast the formulas of the previous section into more 
useful forms. We first compute that (dt/ds)^ = i^i^ + ^l2^ and substitute into (0). Next, 
a simple calculation gives 

^3 = -{ip + cos 6*0) and Vt^ = -(0 + cosO'tp). (13) 

Next, note that t = E^ = g^iCi = sin6'(sin<;/)ex + cos^Cy) + cos6'ez. Explicit evaluation of 
the local Writhe density (^ then gives with (^,p!3D that 

Lk=-— [{^ + ^)ds = — [ ^ \ ds. (14) 
27r J ^' 271 J I + cose ^ ' 

With this last expression, the energy functional (^ is explicitly given in terms of an 
element of the rotation group and its derivatives, as expressed by the angular frequencies 
VL=a and VL^. 

We close this section with a mathematical fine point, which will not affect our 
calculation. Strictly speaking, our configuration space is only locally the group manifold 
50(3). We will exclude the points 9 = it where (0) is singular. Moreover, we need to 
"unwrap" the remaining space. The physical origin of this step is simply the fact that 
rotating the rod by 27r does not return it to an equivalent state, but rather introduces an 
extra unit of Link. Mathematically we simply remember that (p + ip is not to be identified 
modulo 271 (see (p!^), though cj) — ip is. 

V Calculation 

V.A The Path Integral 

We wish to compute the average extension (z) and relative excess Link (Lk) for a twist- 
storing polymer subject to a given tension and torque. To find these properties, we 
must first compute the partition function. At each point along the arclength of the 
polymer, the local orientation will be given by some rotation g. To calculate the weight 
of any configuration entering into the partition function, we simply apply the appropriate 
Boltzmann factor. In the last section we described how the terms of the energy functional 



10 



appearing in this factor can be written in terms of rotations. Using these expressions, it 
is now possible to write down a path integral on the group space: 

1 



Z = / [dg(s)] exp 



(Ebend + ^twist) + 27rr ■Lk+ f ■ z 



(15) 



This partition function gives us the quantities of interest, namely the average chain 
extension (z) and the average excess Link resulting from an applied tension and torque: 



d_ 

Of 



InZ, 



{Lk) 



l_d_ 



InZ. 



(16) 



A direct evaluation of the partition sum in ( |15|) is difficult; fortunately, such an evaluation 
proves to be unnecessary. In this paper we instead extend a standard polymer physics 
trick.B' EH It turns out that the partition sum is closely related to the "propagator" for 
the probability distribution for the polymer's orientation g. We define the unnormalized 
propagator by 

E[g{s)]\ 



^(gf,Sf;gi, Si) 



(sf )=gf 

;{si)=gi 



[dg(s)] exp 



(17) 



The probability Ps{g) for the polymer to have orientation g at position s is then given 
by a multiplicative constant times / dgi \I^(g, s; gi, 0)Ps=o(gi)- More interestingly from 
our perspective, for a long chain log \l/(g, L; gi, 0) becomes independent of g and gi. In 
fact the propagator is then just a constant times the partition function Z. The utility 
of studying the seemingly complicated \I' instead of Z comes from the realization that 
obeys a simple differential equation. We will derive this equation in section |V.B| . Its 
solution for large L is dominated by a single eigenfunction of the differential operator. 
Armed with this knowledge, we will compute in section |V.(J| quantities such as the average 
extension (z) and linking number (Lk) by substituting for Z in the thermodynamic 
relations ([I6|). 

V.B The Schrodinger-Like Equation 

The next step, then, is to determine the differential equation obeyed by \I'(g, s; gi, 0) as 



a function of s. To do this 
length e: 



51 



vE'(gf, Si + e; gi, 0) = J dg^ ^ J^^^^^^^^^^ [dh{s)] exp 



consider the evolution over a short backbone segment of 

5E[h{s)]- 



li(sf )=gl 



k^T 



^(gi,Sf;gi,0). (18) 



Here 5i?[h(s)] is the elastic energy of the short segment of rod from Sf to Sf + e. We 
introduced a normalizing factor to get a continuum limit: as long as this factor does 
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not depend on / or f it will not enter the quantities of interest (see (|T6])). In this 
subsection we will compute the functional integral in ([18D , retaining terms up to first 
order in e, and hence compute d\l//dsf. 

As e 0, we will see that only matrices gi close to gf produce appreciable 
contributions to the path integral. It is therefore possible to write gi uniquely in the 
form gi = exp{—Ta\-a)E,f- Moreover, over the short segment under consideration we may 
take h(s) to interpolate between gf and gi in the simplest way: 

his) = exp ( '~''~^ T^L^ gf. (19) 



e 

The functional integral then reduces to an ordinary integral over T: 

/k...^. tdh(.)] - / exp d3f . (20) 

We have suppressed an overall constant, absorbing it into in (|TBp. The exponential 
factor on the right side gives the invariant volume element of group spaceS near the 
point gf . In the end, this factor will not modify the differential equation that we develop, 
but it is included here for completeness. 

The energy functional 5i?[h(s)] can now be evaluated on the arclength slice of length 
e. With the useful abbreviation 

Maigi) = (gfLggf^^U) = (sin6'sin?/',-sin6'cos?/',cos6'), (21) 

we get that Qa = —Ta/^ and Cl^ = —M-T/e which are constants (independent of s) over 
the short segment. Thus 

5E\h(s)] A,^o C^o (r^ TiMi+T2M2\ ^ ^ 

The factor q-^^/i^bT -^gigj^^g Qg,c\i path from gi = exp(— T ■ Z)gf to gf; as e ^ it indeed 



kills all those gi which wander too far from gf, i.e. all deformations where Ta>^e/ A. 

We also need to express ^I/(gi) in terms of T. Here and below we abbreviate 
\E'(gf, Sf, gi, 0) by ^E'(gf). Define the left-acting (body-fixed) derivatives J'a via 



9^ 



(23) 



dgpi 

and similarly the right-acting (space-fixed) derivatives J7i. Then ^(gi) = e~"^"'^"\E'(gf) or 

^(gi) = ^-f-J^ + ^T^TpJ^Jp * + ■■■, (24) 
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where we abbreviated still further by omitting the basepoint gf on the right-hand side. 

We can now combine ([T8|j20|j2^j2^) and perform the Gaussian integral d^T. 
First complete the square, defining T3 = wC/2e(T3 + ef/C) and = 



A/2e {Ta + efMa/A{l + cos6)), a = 1,2. Choose the normalization N so that the 
limit e — s> reproduces Collecting all order-e terms and using Mi^ + = sin^ 9 
then gives 



COS^^' 



1 1 1 

2" I C ^ 11+ cos9 



f cos9 



T 



A{1 + COS 



1 



C A 



1 



2 \A 



C 




(Ml Ji + M2J2 + M3J3 + J3) 



(25) 



Further consolidation then gives \E' = — (Ti + Sq)"^, where 

and the differential operator 7i is defined by 

1 




(26) 



n 



K 
~A 



2K 



f2 



cos( 



cosy — 



r 



A _ In 1 ^ 



AK 1 + cos 61 
f 1 — cos ^ 



-(- 



4:K1 + COS 6^ 



( J3 + J., 



(27) 



We have arranged the terms in (^) to facilitate a systematic expansion in powers of K^^, 
where = \l Af - f^/A. 

An important property of 7i is that it commutes with both the operators jTs and J^. 
The physical meaning of this property is simply that a uniform rotation of the rod about 
the constant axis changes nothing, and (by the rod's isotropy) neither does uniform 
rotation of the rod about its own t-axis. 

Thus the unnormalized propagator ^ obeys a differential equation which is of 
Schrodinger type, in imaginary time. The derivatives Ja correspond to i/h times the 
usual angular momentum operators, and so on. In the next section, we will exploit the 
quantum mechanical analogy to find solutions to this equation which will in turn allow 
us to determine the quantities {z) and (Lie). 



V.C Solution and Results 

It is now possible to make a direct connection between the eigenvalue problem associated 
to (pTf ) and our polymer problem. 
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In ordinary quantum mechanics, the solution to the Schrodinger equation for a 
symmetric top can be written as a superposition of Wigner functions:l22 



jmk 



-if. 



jmk 



t 



mk vSJ 



(2^ 



Here m and k are angular momenta associated with the operators and J7z, and 
£jmk is the eigenvalue associated with the Wigner function V-!^^. The coefficients Cj^k 
characterize the initial state at time t = 0. 

It may seem difficult to apply ( PS]) to our statistical problem, since in our case jTs 
and J7z are real, antisymmetric operators with no basis of real eigenvectors. Similarly, 
and unlike the case of the wormlike chain, Ti has no particular symmetry. A little 
thought shows, however, that these are surmountable problems. Since one end of our 
rod is clamped, the initial probability distribution \E'(g, 0) may be taken to be a delta- 
function concentrated on g = 1, the identity matrix 9 = ip + cf) = Q. This \E' is indeed 
an eigenstate of jTs — with eigenvalue m — n = 0. The other end of the rod may 
also be considered clamped to 6' = 0, but since we work in the fixed-torque ensemble the 
overall rotation ip + (p is free to take any value. In other words, after evolving \E'(g, 0) 
to \E'(g, L) = e~^^°'^^^'^{g,0) we need to project it to the eigenspace with J^s + J^z = 0. 
Since as noted earlier jTs and jTz both commute with 7i, we may perform the projection 
on ^(g, 0) instead. 

Thus for our problem we should simplify (^) by setting jTa = and J'z = 0, obtaining 
the differential equation that appeared in earlier work:E3 H 4' = — (7Y + Sq)"^, where 




cos ( 



(29) 



K = ^jAf- f 2/4 , (30) 

and So = —{f+f^/2C). The major difference between this equation and that obtained for 
ordinary (non-twist storing) polymers is that the long-wavelength cutoff is now controlled 
by K instead of \J~Aj . 

The operator in (|2^) really is symmetric, and hence will have real eigenvectors (modulo 
a subtlety discussed in appendix B). The solutions to our Schrodinger-like equation will 
then have the form (^) with it replaced by arclength s. For a sufficiently long chain, the 
lowest "energy" solution will then dominate ^ . The thermodynamic properties of the 
polymer can then be determined by remembering that \E', the unnormalized propagator, 
becomes equal to a constant times the partition function Z, and applying (p!6D . 
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We gain further confidence in the above analysis when we note that the terms set to 
zero in (|27|) include some which are linear in the applied torque r. For reasons outlined 
in section PTT| , we do not expect these terms to play a role in the determination of the 
lowest energy eigenvalue. The model that we defined is non-chiral and therefore cannot 
tell the difference between over- and undertwisting. 

We must now compute the lowest eigenvalue of the differential operator in (PU]). 
Finding it would be a straightforward task were it not for the singularity in the potential 
term when 9 ^ n. This singularity is associated with the backbone tangent i looping 
around to point anti-parallel to the end-to-end displacement vector +ez. Physically, this 
situation corresponds to the onset of supercoiling. When the applied torque is too high 
or the tension is too low, the chain will begin to loop over itself. Since real chains cannot 
pass through themselves, they begin to form plectonemes. In our phantom chain model, 
there is no self-avoidance, and so the chains can pass through themselves, shedding a unit 
of Lk as they do. The mathematical pathology associated with the ^ — > tt singularity in 
(p9D is therefore an inevitable consequence of our model's neglect of self-avoidance. 

The physical breakdown of the phantom chain model and the corresponding 
mathematical problem of the 6' — > vr singularity can be avoided by assuming that the 
backbone tangent t remains nearly parallel to the -t-Cz-axis. Such a situation is indeed 
realistic for a chain under sufficient tension, or more precisely, for a sufficiently large K 
(POI). In this regime, we can then perform a perturbative expansion about 6 = 0. The 
singularity of ( pQ]) does not affect low orders of perturbation theory. The singularity 
can still enter nonperturbatively via "tunneling" processes, in which the straight 6^0 
configuration hops over the potential barrier in (PU]), but these will be exponentially 
suppressed if the barrier is sufficiently high, a condition made more precise in appendix B. 
The perturbative regime is experimentally accessible: we will argue that it corresponds 
to the solid symbols on Figure ||. Outside this regime, the phantom chain model is 
physically inappropriate, as explained above, and so a full nonperturbative solution of 
our model would not be meaningful. 

We can simplify the problem by changing variables from 9 io = 2(1 — cos 6'). In 
terms of p the spherical Laplacian J"^ = -^^^ sm9-^ becomes (1 — p'^/A)dp^ + (1 — 
3pV4)p-i9p, so 



1^ 




p2 f2^4 



16 - 4p2 



(31) 



where = p~^dppdp. We have not made any approximation yet. 

We now construct a perturbative solution to the eigenvalue problem defined by 
(|3lD. In the quantum mechanical analogy this equation describes a two-dimensional 
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anharmonic oscillator, with p interpreted as a radial coordinate; thus the problem can 
be solved using the method of raising and lowering operators. 
Switching to Cartesian coordinates, we set 



Now (^) can be rewritten as 7-^ = 7-^o + ^'^i where 
Ho = ^(A^ + ATb + l), and 



5n 



K 
~A 



(33) 



m\ 4 

Here Na = A+A- and Mb = correspond to the usual occupation number operators 

in the quantum mechanical analogy. It is now straightforward to calculate the lowest 
energy eigenvalue as an expansion in to obtaint 
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Remarkably, this is exactly the same formula as the one appearing in the wormlike chain 
model; the only difference is that K is now defined by (|30|) instead of by \J Af. The 
last two terms retained will now give anharmonic corrections to the simple lowest-order 



1 c _, 

calculation announced earlier The ellipsis represents terms of higher order in K~ than 
the ones kept. We explore the status of such terms in appendix B. In particular, the last 
term of (pT)) has been dropped altogether. Since the expectation value of this term is 
obviously divergent at p = 4 {i.e. the antipode i = —z), a certain amount of justification 
will be needed for dropping it. 

From this eigenvalue, the mean extension and the average linking number for a given 
tension and torque can be found using (0), (^), and (|3^): 




More accurate versions of these formulae are given in appendix B. By solving the second 
of these equations for the torque we obtain the new, effective twist rigidity Ces by noting 
that f{f,Lk) « {27rLk/L)C^s{f) + 0{K-^), where Ces{f) is given by the formula ®. 
This formula describes the "thermal softening" of the twist rigidity alluded to earlier. 
The effective rigidity Ccs{f) is reduced from the bare, microscopic value by a factor which 
arises from thermal fluctuations. 
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Combining (|35|) and (^) together with the definition of K in (|30D produces a formula 
for the average end-to-end extension for a polymer subject to a linking number constraint 
and an applied tension. In section |VI| we will compare this theoretical prediction to the 
experimental results of Strick et a/0 and AUemand and Croquette.l^ 

V.D Onset of Non-Perturbative Corrections 

The theory described above is only valid in the regime where the phantom chain model is 
appropriate. In this section, we extend our analysis by estimating the effect of plectoneme 
formation close to its onset. As discussed above, our model is unable to include 
such effects quantitatively, as it lacks the self-avoidance interaction which stabilizes 
plectonemes. Instead, in this section we will suppose that the main consequence of the 
singularity is to allow each segment of the polymer to be in one of two configurations. 
In the first instance, the polymer fiuctuates about a nearly straight conformation and 
can therefore be described by the theory developed in the preceding sections. In the 
second instance, the polymer is driven across the "tunneling" barrier into a standard 
kink conformation as depicted in Figure 0, gaining approximately one unit of Writhe. 
Our improved formulae will have no new fitting parameters beyond the ones already 
introduced. 

We are interested only in the initial stages of plectoneme formation and so it will be 
sufficient to approximate each plectonemic coil by a circle. The energy required to form 
such a loop is 

Here the first two terms represent the energy costs associated with bending the polymer 
and contracting against the imposed tension. The last term gives the elastic twist energy 
released as Twist gives way to Writhe. Maximizing the energy release, we find the optimal 



radius of a coiled segment to be = J A/2f, so that for f > the presence of a kink 



lowers the energy of the polymer by ^E^/kj^T = 2'K{\j2Af — f). 

We now imagine the polymer to be made up of segments of length 2itR. Each of 
these segments may be in the extended or the plectonemic kink configuration. Actually, 
we will consider two possible types of kink: one in which a unit of Twist is shifted into 
Writhe, and its mirror image which generates a negative Writhe as well as a counteracting 
positive Twist. The reverse kinks are energetically unfavorable for appreciable applied 
torque, but we retain them to eliminate any asymmetry in the excess linking number. 

For the modest applied torque considered here, it will be sufficient to treat a dilute 
gas of positive and negative kinks. Denoting the population of kinks by n_ and that of 
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reverse kinks by n+, we have 

, , K.L ( AE+\ , , 

Here AE'^. = 2Tx{\j2Af + f) is the energy of a reverse kink, and k is a numerical factor 
of order unity arising ultimately from a functional determinant. Since we do not know 
how to compute k we set it equal to unity. 

The effect of the kink/anti-kink gas is to modify (|35|) and (|36D , producing shifts in 
the average extension and average linking number: 

a{j^ = - (e-^^-/'=B^ + e-^^+/^BT) (39) 



These expressions are to be added to (|35|) and (0). The latter expression can then be 
solved for f to get a corrected version of (H). 

The model proposed here for plectoneme formation is too simplified to give 
quantitative predictions about the non-perturbative regime. However, the model does 
allow us to predict the onset of these effects and confirm that the data we select are not 
affected by plectoneme formation (see Figure |ip. 

VI Fit Strategy and Results 

The extension function {z{f,Lk)) derived in the previous sections describes an achiral 
elastic rod. Before making direct comparisons of this formula to experimental data, we 
will extend the model somewhat. So far we have neglected structural changes in the DNA 
at a microscopic level. In particular, we have omitted effects related to the intrinsic 
stretching along the polymer backbone. Recent experiments have investigated these 



effects;!^^' 1^ in particular, Wang et al. found a small change in the relative extension 
of //7, where 7 = 1100 pN is the intrinsic stretch modulus. For moderate forces we 



may simply add this shift to the extension formula found in the previous section.@ El 
For the highest forces we consider (8.0 pN), this translates into a relative extension of 
about 0.007, which is hardly noticeable in Figure Nevertheless, we will include this 
correction as it improves the quality of our fit slightly without introducing a new fitting 
parameter. 

In addition, we will also consider the possibility of elastic couplings which do not 
respect the inversion symmetry of the model that we consider. In reality the DNA we 
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seek to describe is chiral, and so at some level we expect this fact to show up as an 
asymmetry between overtwisting and undertwisting in Figure ^ One way that chirality 
might enter a model for DNA is through an intrinsic twist-stretch coupling.0 EH 
This coupling results in a change in relative extension of —kjiTDujQa/^, where D is the 
twist-stretch coefficient. The near symmetry in the data of Figure ^ indicates that the 
effects of such a coupling will be small in the region of interest. Although the coefficient 
D will turn out to be comparable in size to the bending coefficient A, the shortening 
due to bend fluctuations dominates that due to the elastic twist-stretch coupling. This 
disparity arises because bend fluctuations are diverging as K 0. 



As mentioned in section |T|, an anisotropy between the "tilt" and "roll" elastic 
constants coupled together with the associated twist-bend coupling term might also 
produce an asymmetry between positive and negative a. We investigate this possibility 
in appendix A and flnd that the corresponding chiral entropic elasticity terms are not 
measurably different from the twist-stretch model over the range of stretching forces 
studied. 

Putting the intrinsic corrections associated with 7 and D together with the 
perturbation theory result of the last section, we obtain a theoretical prediction for 
the relative extension as a function of applied force and overtwisting. For the purposes 
of comparison to experiment, we will now switch from the variable Lk to the relative 
overtwist a which is deflned with respect to the helical pitch of DNA: a = 27TLk/ujQL. 
Then, 




L / \ ^ kj,T 4 32 / \L/ 7 



Formula ( ^T]) is our flnal result for the high-force (or more precisely, large K) extension of 
a twist-storing polymer subject to a torsional constraint. Here A{z/L) is the expression 
in 



). To compare our result to the experimental data,ll3 Il3 we solved (P^, and 
for f in terms of /, cr, then substituted f and K into (^). 
Apart from the intrinsic stretch and twist-stretch terms described above, (|4l|) contains 
two additional small reflnements. One of these appears in the last term, where flnite-size 
effects have been accounted for. This term can be understood by writing the extension as 
an expansion in terms of the transverse components of the backbone tangent. Deflning 
the complex variable a{s) = i{s) ■ (e^ + iCy) and its Fourier components Op, we have to 
lowest order 

l-^E(l«/) + --- • (42) 



L/ 2 p 
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The leading entropic reduction of {z) in (|35D is then easy to evaluate, including finite- 
length effects. As the main effect of an applied torque is to decrease the effective force 
and change the low wavenumber cutoff in our theory from \J Af to K, we know how to 
modify the usual tangent-tangent correlation function to yield 

4A °° 1 

E(l«/) = 



This expression should be compared with the leading-order correction obtained from the 
infinite-rod calculation in section 0: 

i2\ 2A 1 

« ^- (44) 

The difference between the two terms is 2A/LK^. To obtain the finite length formula we 
must subtract this difference from the result obtained in the last section; the resulting 
correction appears in the last term of (^Tp. Note that for the restricted values of K that 
we consider (see below), this contribution to z/L never exceeds 0.002 for the data set we 
analyze. 

The other refinement introduced in (14T|) is that for convenience we replaced 
{1/2K){1 + 1/64:K^) by {K^ - l/32)~^/^ Since we will restrict our fit to K"^ > 3, the 
difference between these expressions is negligible. Finally, in Appendix B we give even 
more elaborate versions of (|36[) and (|4l|), in which higher-order terms of perturbation 
theory have been retained; these corrections are small though not negligible at low forces. 

We have now established an expression for the mean extension as a function of applied 
tension and torque. Using the ENS group's dataP m we fit this formula (actually, the 
more accurate one given in appendix B) to determine the parameters in our model: the 
microscopic bend persistence length A, twist persistence length C, twist-stretch coupling 
D, and polymer arclength L. Of these parameters, only C and D are really unknown; 
A has already been measured in other experiments, and L can be determined from the 
points with a = using the ordinary worm-like chain model. The agreement between 
our best fit value of A and earlier experiment£3' 0' HI serves as a check on the theory. 
Other parameters appearing in (|4l|) , namely Uq = 1.85 nm~^ and 7 = 1100 pN,0 are 
independently known and are not fit. 

The least squares fit was performed using a gradient descent algorithmlH in the 
parameter space defined by A, C, D, and L. The best fit was obtained for A = 49 nm. 
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C = 109 nm, D = Q7 nm, and L = 15.6/xm. Here L is the length of the construct from 
AUemand and Croquette's experiment.0 The corresponding length for the Strick et al. 
data setlH was determined separately using the a = points from that set and was not 
fit. In all, 69 data points from the experiments of Strick et a/." and of Allemand and 
Croquette^ were used in the procedure. The data points were selected based on three 
criteria. The first cuts were made on physical grounds. It is known that for high applied 
forces (/ > 0.4 pN) DNA undergoes structural transformation or strand separation when 
a < —0.01 or (T > 0.03 (D. Bensimon, private communication); here of course we cannot 
use linear elasticity theory. We therefore omitted such points from the right side of 
Figure |l|. (No points were omitted from the left side.) To avoid biasing the data, in 
the fit we excluded the symmetric region |cr| > 0.01 from the set of points used with 
/ > 0.4 pN. 

The second set of cuts was applied for mathematical reasons. Our perturbative 
expansion is in powers of K~^: we required > 3. We discuss this choice in appendix B; 
for now we note that perturbation theory produces excellent agreement with experiment 
for the wormlike chainO even for K > 1. Choosing > 3 eliminates all of the / = 0.1 
and 0.2 pN data points from our fit. To confirm that we were being selective enough, we 
tried other values of the threshold (between 2.5 and 4.5). This action did not significantly 
alter our fit results: in every case we found C > 100 nm. 

Finally, in addition to these two sets of data cuts, we also imposed a "tunneling" 
criterion described in appendix B: the idea is to ensure that the lowest energy eigenvalue 
of the operator Tio in ( |33D is smaller than the barrier that restrains the system from 
falling into the unphysical singularity. 

The reasonable agreement in Figure |I| between our theoretical curves and the data 
outside the region we fit (including the 0.1 and 0.2 pN curves) indicates that our choice 
of cuts is a conservative one. As a further check, the dashed lines in Figure |l| show our 
fitting function without the non-perturbative correction described in section |V.D| : we see 
that these lines do not deviate from the solid lines in the range of data we retained. 



VII Discussion 



The global fit shown in Figure |l] indeed resembles the experimental data. The least 
squares fit determined the bending stiffness, the twist rigidity and the intrinsic twist- 
stretch coefficient of DNA. As stated earlier, the fit to the bending rigidity produced the 
known value and thus serves as a check on the theory. The chiral asymmetry is a small 
effect, and so the available data do not afford a precise determination of the twist-stretch 
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coupling D. Thus our fit is mainly a measurement of C. 

The twist rigidity obtained by the fitting procedure is somewhat higher than what 
earlier experiments have found (see section [TTTD. We cannot give a quantitative estimate 
of our fit parameter errors, since some of the datali^ do not have error bars, but we 
note that forcing C = 85 nm or less gives a visibly bad fit. One might worry that this 
discrepancy was due to some sort of failure of perturbation theory, despite our great 
care on this point. The fact that we keep finding large C as we tighten the data cuts 
gives us additional confidence on this point. Similarly, our large value is not an artifact 
of DNA denaturation induced by tension, since that would lead to a spuriously low fit 
value. There remains the intriguing possibility that on the contrary, imposed tension 
suppresses spontaneous local denaturation, increasing the integrity of the DNA duplex 
(J. M. Schurr, private communication); in this case our large C more accurately reflects 
the linear elasticity than the other, lower, values. 

The discrepancy with earlier work may be more apparent than real, however: if we 
do not allow for a tension-dependent thermal reduction of the twist rigidity as in (§) and 
instead fit the data to a constant twist rigidity, then we obtain C^s = 82 nm, a value 
closer to those found in the other experiments.^ The quality of this fit, obtained with a 
tension- independent rigidity, is slightly poorer. In any case, a large value of C/A is not 
paradoxical and in particular need not imply a negative Poisson ratio for our model's 
rod: random natural bends in DNA reduce the effective bend stiffness A measured in 



stretching experiments, but not Cp^ and so the ratio of C to the true elastic bend 



stiffness is closer to unity than it appears from our effective-homopolymer model. 

y I9ni 

Recently, Bouchiat and Mezard^-^l have also determined the twist rigidity of DNA 
using the experimental results of Strick et a/.Hll They derived formulae equivalent to ( |5BD 
and (|36|). Then using an exact ground state solution to a cut-off version of (0), they 
reproduced the observed extension curve {z{f,a)) in Figure ^ over a wider range than 
we have shown. The result of this calculation is a ratio of C/A of approximately 1.7. 

While both approaches are similar, our perturbative approach precludes us from 
analyzing the lowest force curves that Bouchiat and Mezard discussed. As described 
above, we excluded these data because we expect physical difficulties with the phantom 
chain model in this regime; the same difficulties, it would seem, apply to the analytical 
results of Bouchiat and Mezard. In particular, at small applied tension, the backbone's 
tangent vector t will wander from the z-axis. If it wanders too far, the system will be able 
to see through the tunneling barrier to the singularity; or in other words, the results will 
be corrupted by the failure of Fuller's formula for Wr. Bouchiat and Mezard approached 
this problem by introducing a new intermediate-length cutoff 6 = 6 nm into the problem. 
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The physical meaning of this cutoff in terms of the mechanical properties of DNA is not 
clear to us. Moreover, taking it to be 2.5 nm or less spoiled the simultaneous fit at all 
values of /. 

In contrast, our perturbative treatment avoids the singular-potential problem 
altogether by restricting to a regime where the phantom chain model is valid. Our 
model has no extra scale corresponding to b, and yet fits all fixed-force curves in its 
domain, in two different experiments, with one value of C. 

In their paper, Bouchiat and Mezard also gave Monte Carlo results. Earlier work by 
Marko and Vologodskii has also taken this approach.0 Here it is possible to implement 
self-avoidance, though knot rejection is still difficult. The advantage of analytic formulae 
such as (^) is that they permit global, systematic least-squares fitting of {z{f,a)) to 
the data. Moreover, for practical reasons Monte Carlo simulations must again impose 
a short-distance cutoff of at least several times the DNA radius, unlike our analytical 
approach. 

VIII Conclusion 

In this paper we have investigated the statistical mechanics of a twist-storing polymer. 
This type of molecule differs from a traditional polymer in being unable to relax out an 
applied excess Link. When such a chain is left unconstrained, the twist simply decouples 
from the bend fluctuations. The thermally accessible conformations are then identical to 
those for an ordinary polymer. In the case that such a polymer is subject to a torsional 
constraint, however, there will be a coupling between the bend fluctuations and the twist. 
It is this coupling that we have investigated. One of our goals was to show how single- 
molecule stretching experiments can provide a new window onto the nanometer-scale 
mechanical properties of DNA. 

Due to the complications associated with self- avoidance, we considered only chains 
held nearly straight by tension, then analyzed the statistical mechanics of the resulting 
"torsional directed walk" . We mapped the polymer partition function onto the solution of 
a Schrodinger-type equation for the orientation distribution function. From this solution, 
we were able to find the entropic extension and the overtwisting of a polymer subject to 
a tension / and relative Link excess a. 

The theory we developed quantitatively reproduces the results of supercoiled single- 
molecule DNA stretching experiment (see Figure |1]). The agreement was achieved 
by fitting the twist persistence length, yielding C = 109 nm. The large twist rigidity 
differentiates DNA from traditional polymers and makes possible the coupling of the 
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twist and bend degrees of freedom that plays a central role in our theory. 

Apart from reproducing the experimentally observed physics, our formulae make 
another prediction: the twist rigidity is renormalized (see ([36|)). The effective rigidity 
Cesif) is a function of the applied tension. According to (§), it is hardest to twist the 
polymer when it is pulled straight; this is the bare, microscopic stiffness. It is the same 
rigidity that resists twist at the shortest length scales, and so enters the energetics of 
structures such as the nucleosome. As the tension is relaxed, thermal fluctuations begin 
to play a role. Now when a torque is applied, the polymer does not resist as much; the 
bend fluctuations have softened the torsional rigidity by absorbing some of the imposed 
excess Link. As discussed above, this phenomenon is purely thermal; no such effect 
appears in the linear elasticity of a macroscopic beam for small applied torque. 

If one naively extends this thermal effect to zero tension, one sees that the torsional 
rigidity vanishes completely. Of course, our phantom chain model precludes us from 
considering this case; however, other recent wor k0 has considered this related problem 
using an explicit self-avoidance term: indeed the effects of a torsional rigidity do 
become unimportant to the behavior of twist-storing polymers at zero applied tension 
or, equivalently, at extremely long length scales. 



Appendix A: Chiral entropic elasticity 

In this appendix we introduce an additional element of realism into our model, namely 

the intrinsic helical pitch 27t/ujq. For DNA this pitch corresponds to ujq = 1.85/nm. The 

helical structure breaks the inversion symmetry of the problem by allowing two additional 

I97I 

terms in the energy functional.!^ In principle these explicitly chiral terms could introduce 
an asymmetry between overtwist and undertwist into our results. We will find this chiral 
entropic elasticity and show that it has a different dependence on stretching force from 
the intrinsic twist-stretch effect discussed in section Thus in principle the two effects 
could be distinguished experimentally. 

In this appendix we are interested in chiral effects, manifested by odd powers of a in 
the extension z{f,a), in a model of DNA without intrinsic stretching. We will see that 
such terms are small. Hence we can use a simpler calculation than the one in the main 
text: we will drop 0(o"^) and higher, and we will use the Gaussian (or equipartition) 
approximation to the statistical sums. Since odd-power terms are completely absent in 
the achiral model of section |V| above, we can simply add the ones we find to the results of 
that model to get a leading approximation to the full chiral entropic elasticity formula. 

Another approximation we will make will be to drop terms suppressed by powers of 
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l/cuo, since this length scale is much shorter than both the persistence lengths and the 
scale \J A/ f of important fluctuations. 

As discussed in section |TTT|, chirality can enter through the anisotropic bending 
rigidities associated with the "roll" and "tilt" axes of DNA monomers. In this appendix 
we will choose a material frame different from the one in the main text: here our frame 
rotates with the intrinsic helical twist. This choice is convenient in that the anisotropic 
elasticity appears constant in this frame: (|I]) becomes simply 

F 1 i-L W ^ rL 



k^T 2 



/ dsUn^ + ^,^^), and ^ = U C'Qs'ds. (45) 

Jo ^ ' Kjil 2 Jo 



Here we have introduced two microscopic bending constants, A[ and Ag. Now even 
the unstressed state will be chiral: as the body-fixed frame {Ei, E2, E3 = i} rotates 
uniformly at frequency ujq along the polymer, it turns the bend anisotropy with it. Since 
the £'i-axis corresponds to the short axis of a basepair, we expect A[ > A'2. 

Apart from the bending anisotropy, the symmetries of DNA admit an explicitly chiral 
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term associated to a twist-bend coupling with coefficient G.lfil With these two terms, the 
mechanical-equilibrium state of the stressed molecule will no longer be given by the 
uniformly twisted configuration. Instead, we make an ansatz for a new helical ground 
state: go = exp(^Li) exp(c<JsL3), to be justified below. Here u includes a finite piece 
associated with the rotation of the unstressed molecule, so that uj = ljq{1 + a). The 
small angle ( remains to be determined by the condition of mechanical equilibrium. The 
elastic energy functional for the model is then given by: 

E 1 fL 



ds {A[n,^ + A'^n^^ + C\n^ - uo)^ + 20^2(^3 - ^0)} - f-z. (46) 



In contrast to the discussion in the main text, in this appendix we will work in the fixed- 
Li ensemble. Thus we do not need any Lagrange multiplier associated with the Link 
constraint. 

It will prove convenient to introduce the combinations A = {A[ + ^2)/2 and 
A = {A[ — A2)/2. We emphasize that A is not necessarily equal to A from the coarse- 
grained model (|l|); the exact relationship will emerge in due course below. The chiral 
terms that couple to the intrinsic helical frequency Uq are then proportional to A and G. 
Note that i > 0. 

We can now determine the helix angle ( characterizing the mechanical-equilibrium 
state. First write a small fiuctuation from go as g(s) = g(s)go(s) with go as above and 
g(s) = Q-'^o'i^)^'^ , Substituting into (|ll],|l2|) then yields the f2.j's. Setting the first variation 



of (|46|) to zero then yields three equations expressing the condition that go be the stressed 
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mechanical-equilibrium state. One of these selects (: 

C = = ~ — rr- (47) 

The other two are satisfied trivially, justifying our ansatz for go- In deriving the above 
relations we used the fact that we are working in the fixed- cr ensemble. Thus the boundary 

— * 

conditions clamp the rod at both ends, fixing T — there, and so we may discard total 
derivative terms. 

For illustration, and to keep the calculation simple, we will now make the additional 
assumption that the chiral parameters A, G are both smaller than A, C, and accordingly 
work to leading nontrivial order in the former. Wc can then easily diagonalize the part of 
the energy involving the latter using Fourier modes. Setting (Ti(s)-|-iT2(s))e''^* = J2&'^^aq 
and T3(s) = Ee^'^'^g (note that = 0*) yields 

^ -/•L + eo + 6 + 6, (48) 



where 



^0 = ^^[{Ap'-u;{C'a + 2G0p + f)\a,\' + C'p'\cl^, 



2 p 



= ^ [^(u; -p)p + C (A[u;{uj - p) + f - C'ujp]] {(ppa^^-p - c.c.) 



2 p 



-p{p - 2uj) —p 



{oLpOLlui-V + c-c.) . (49) 



In the above formulae, a; = a;o(l + cr) gives the angular frequency for the stressed minimal- 
energy state. The sums are for —oo < p < oo (the physical short-scale cutoff will prove 
immaterial). As mentioned above, we will treat ,^12 as perturbations to ,^o- 
In the harmonic approximation, the mean extension has the simple form 

^ %lnZ^l-WM + ---- (50) 



L/ Ldf 2^ 

We define D{p) = L{\ap\'^) and compute this two-point correlator perturbatively. 

The unperturbed Dq{p) is obtained via equipartition, or equivalently by performing 
the Gaussian (harmonic approximation) functional integral over ctp and a* in ^0, yielding 

A{p^ - 2qop) + / 
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where 

?o=(C"-^)^. (52) 




The next step in determining D{p) is to calculate the first two corrections, Ili{p) and 
112 (p), induced by ,^1 and ,^2 respectively. We define these as lowest-order corrections to 
the full two-point function: D{p) = Dq{p)[1 + Dq{p)(JIi{p) + Il2{pj)]. Start by expanding 
the contribution e~^^ to the Boltzmann factor. There is no first order correction, so we 
go to second order: 



{A{2uj - p) + GioCf 



U2(p) = 

A{i2co - py - 2qoi2u - p) + f 
^ 2P + CAG 

~ P J • (53) 

The second correction arises from the expansion of the energy in powers of ^1. Once 
again we go to second order: 

As mentioned above, we have dropped terms of order o"^ and higher: only odd-power 
terms will create chiral corrections to the extension curve, and we content ourselves with 
investigating the linear ones only. Putting the results of (Q) and ( |53D together gives the 
propagator 

D{p) = [Doipy - U,{p) - Mp)]'' . (55) 

To get (^) we summed chains of Gaussian graphs, similarly to the random-phase 
approximation in many-body theory. 

The relative extension can now be computed from (|50|): 



1 

1-— / dpDip) 



= 1-1(A/(1 + Fa))"'^'. (56) 

In this formula we have identified A = A~2A^ / A — G"^ /C as the effective bend constant, 
coarse-grained over a helix turn. (Had we kept O(cr^) terms we could have made a similar 
identification of the coarse-grained twist constant G in terms of A, A,G' ,G.) We also 
defined 

9n^ (A A'.\ 

(57) 
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The key observation is now simply that F in (|57D is positive. 

Thus we have found a chiral entropic elasticity effect: the formula of the main text 
for 1 — z/L get muhiphed by the asymmetric correction factor {1 — Fa/2). (This analysis 
corrects an erroneous claim^ that no such factor exists.) 

The dependence of this chiral contribution to z/L on the stretching tension is different 
from the intrinsic twist-stretch term introduced in the main text, equation (PT|), and so 
in principle the two effects could be disentangled by fitting to data. In practice, however, 
the chiral effect in Figure |I| is too small to make any definite statement. Instead we tried 
eliminating the D term in (^) and replacing it by the F term in (|5BD, which yields an 
equally good fit but with F = —1.6. Since this value is not positive, contrary to the 
prediction in (^), we conclude that the twist-stretch coupling D is needed to explain 
the asymmetry of the experimental data. This conclusion is qualitatively consistent with 
an earlier analysis!! of the highest-force data; here the chiral entropic effect is very 
small (see (|56D ). Encouragingly, the fit values of A,C are similar to those quoted in the 
main text — our measurement of C is not sensitive to the precise mechanism of chiral 
symmetry breaking. 



Appendix B: Domain of validity 

In this appendix we endeavor to justify our perturbative approach to torsional directed 
walks, and in particular establish its domain of validity and hence the subset of the 
experimental data which falls into that domain. 

Tunneling 

As we have mentioned several times, the Schrodinger-type equation defined by (|29D suffers 
from a singularity at the antipode 9 = n. Indeed, the operator H has no eigenstates at 
all. We have emphasized that this singularity is caused by our unphysical omission of 
self- avoidance effects, but it is still necessary to have some criterion for when the details of 
the nonlocal interaction correcting the problem will be unimportant, and some practical 
scheme for calculating in this regime.^ ^ ^ 

The key point to note is that if we let t = and imagine solving our problem for 
negative (unphysical) values of t, then our problem disappears. Analytically continuing 
the ground-state eigenvalue in the complex t-plane back to positive (physical) t yields 
a result which is finite but no longer real: for small t its imaginary part gives the 
probability of a rare barrier penetration process. The real part is an approximate 
eigenvalue describing the metastable state and controlling the intermediate asymptotics 
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of this is the number we seek. When the imaginary part is small, the real part can be 
obtained from the lowest orders of perturbation theory, even though eventually at high 
orders the series diverges. 

We can estimate the imaginary part of the eigenvalue by finding the saddle point (or 
"instanton" or "bounce" or "domain wall" solution) of the functional integral giving 
rise to (^). This is the function 6{s) satisfying the ordinary differential equation 
A^e = dV/de, where V{e) = {1 ~ cos 6) {K^ - t{l - cos 6) / {1 + cos 6)) . The elastic 
energy of this configuration is then given by 



E 



ds 



1 



deV2V, 



(58) 



where 6i is the "turning point" , where V^(^^i) = 0. The imaginary part of the analytically- 
continued eigenvalue is then proportional to , Numerical evaluation shows that 
this factor is smaller than 0.02 when t < 0.6(^4/ — 1.6), and we have imposed this as one 
of the conditions selecting the data points used in Figure ^. 



Perturbation theory 

From the previous subsection and the references cited there we know that when the 
tunneling criterion is satisfied perturbation theory will be an asymptotic expansion, which 
we may approximate by its first terms. In this subsection we will quote the eigenvalues 
of (|31| ) obtained using second-order perturbation theory. In the last term we expand 
p''/(l — p^/4) in power series, since each succeeding term is formally suppressed by a 
power of K^^] we keep the terms + p^/4:. We again abbreviate t = 

Using the operator notation of the main text, we find H = Hq + 6TC, where Hq is the 
first line of (|3^) and 



sn = — 



1 

8A 



1 / t 



(A' + 3A+^B+^ + 3A, 
9t 



2 9 



)(A' + 2A'^4 
2(1- 



t 3t 



(59) 



plus terms annihilating the perturbative ground state. From this we compute zeroth 
through second-order shift: 



S 

Wi 



t 



Wo 



1 - 



AK QAK^ 
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8K 32K^ 
9 ^ 

4^ ^2K^9K2' 



(60) 



Taking thermodynamic derivatives as in the text (see (|16|)), and recaUing t = t'^/A, gives 



z/L 



1 

C 



1 



AAK 



1 21 
2K ^ 



f2 



(2 



15 9 . 
8^ ^ 8^' 



+ M 



—d ^ 



+ 



f2 



16A" 



8^2 



+ M 



where 



M 



f2 



9x/ 15 112 X 



(61) 



(62) 



The corrections for kinks, (|39D and (0), and the other corrections in (|41|) must be added 
to the expressions (|6ll) . The resulting formulae are the ones actually used in the fit shown 
in Figure |I]. 

We are now in a position to state the conditions for perturbation theory to be useful. 
Our expansion is in powers of K'^ and f^/16-ft'^, so both of these must be small. To be 
more precise, we imagine holding the force / fixed while varying the torque f, as in the 
experiment. The coefficient of in z/L then gives the information we need to obtain 
the twist stiffness. Comparing the highest-order term of this coefficient retained above 
to the leading term, we find their ratio to be less than 10% when > 3. This explains 
another of the cuts made on the data in the text. We should also require that f^/16-ft"^ 
be small, but this is automatically satisfied when the other imposed conditions are. 
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Figure 1: Relative extension of A-DNA versus applied force / and overtwist a: a single 
global fit to two experiments. Fitting our model to the solid points shown correctly 
predicts many of the open symbols shown, even though they were not used in the fit. 
On the left are experimental data from Allemand and Croquette:^ from top to bottom, 
the curves are at fixed force 0.388, 0.328, 0.197, and 0.116 pN. The error bars reflect 
the measurement of extension; estimated errors in the determination of the force are not 
shown. On the right are data from Stride et al. from top to bottom, the curves are 
at fixed force 8.0, 1.3, 0.8, 0.6, 0.3, and 0.1 pN (error estimates not available). Points 
corresponding to /, a where the DNA is known to denature or undergo structural change 
have been omitted from the right hand graph. Solid symbols are within the range of 
validity of our model (for example, all solid symbols have > 3, see text); open 
symbols were not included in the fit. A total of 69 experimental data points were used 
in the fitting procedure. Some of these points are not shown; they had force not equal 
to one of the ten values listed above. The solid lines are a single global fit to both 
datasets using the theory developed in the text (see (^T])). The dashed (higher) lines 
are the same theoretical curves but without our estimated non-perturbative contribution 
(section |V.D|) . 
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Figure 2: Diagram showing the ideahzed circular loop model of a plectoneme. The 
twisted and slightly writhed conformation above is shortened by the coil circumference 
as the plectoneme forms. 
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